Code
# Ci-dessous les librairies R utilisées pour cette analyse
library(tidyverse) # pour la manipulation et visualisation de données
library(sf) # pour traiter les données spatiales
library(tmap) # pour générer des cartes
library(aws.s3) # Pour accéder aux données volumineuses stockées en ligne
library(geodata) # Pour avoir les frontières administratives
library(curl) # Pour télécharge
library(lubridate) # Pour manipuler des dates
library(httr)
library(mapme.biodiversity)
library(units)
# Create a dataframe of the coordinates
coordinates <- data.frame(
location = c("Beroroha", "Beronono", "Tsivoko", "Makaykely"),
latitude = c("21°40.504’S", "21°21.669’S", "21°17.712’S", "21°28.074’S"),
longitude = c("45°9.571’E", "45°14.885’E", "45°22.732’E", "45°21.896’E")
)
# Function to parse DMS and convert to decimal degrees
ddm_to_decimal <- function(dms){
# Extraire les degrés, minutes et secondes
parts <- regmatches(dms, gregexpr("[0-9.]+", dms))[[1]]
degree <- as.numeric(parts[1])
minute <- as.numeric(parts[2]) / 60
direction <- ifelse(grepl("S|W", dms), -1, 1)
# Converstion au format décimal
decimal <- direction * (degree + minute)
return(decimal)
}
# Apply function to each coordinate
coordinates$latitude_decimal <- sapply(coordinates$latitude, ddm_to_decimal)
coordinates$longitude_decimal <- sapply(coordinates$longitude, ddm_to_decimal)
# Convertir les coordonnées en objet sf
coordinates_sf <- st_as_sf(coordinates,
coords = c("longitude_decimal", "latitude_decimal"),
crs = 4326)
makay_ap <- st_read("data/Makay_wgs84.geojson", quiet = TRUE) %>%
rename(Statut = SOURCETHM) %>%
mutate(Statut = recode(Statut, "zone tampon" = "Zone tampon"))
# Fusion des différentes aires
makay_union <- makay_ap %>%
st_union() %>%
st_sf() %>%
st_make_valid()
buff_size <- 20
# Create a 20km buffer around Makay PA
buffer <- makay_union %>%
st_buffer(dist = set_units(buff_size, km)) %>%
st_as_sf()
makay_ext_buffer <- buffer %>%
st_difference(makay_union)
aoi <- makay_union %>%
bind_rows(makay_ext_buffer) %>%
st_cast("POLYGON") %>%
mutate(name = c("Makay", "Périphérie 20km"))
tc_2000 <- rast("/vsis3/projet-betsaka/diffusion/PA-impact-on-deforestation/mapme/gfw_treecover/Hansen_GFC-2024-v1.12_treecover2000_20S_040E.tif") %>%
crop(buffer) %>%
aggregate(fact = 4, fun = max)
|---------|---------|---------|---------|
=========================================
Code
lossyear <- rast("/vsis3/projet-betsaka/diffusion/PA-impact-on-deforestation/mapme/gfw_lossyear/Hansen_GFC-2024-v1.12_lossyear_20S_040E.tif") %>%
crop(buffer) %>%
aggregate(fact = 4, fun = max)
|---------|---------|---------|---------|
=========================================
Code
loss <- lossyear
values(loss) <- ifelse(values(lossyear) > 0, 1, NA)
loss_recent <- lossyear
values(loss_recent) <- ifelse(values(lossyear) > 21, 1, NA)
tmap_mode("view")
tm_shape(tc_2000) +
tm_raster(title = "Forest cover in 2000 (%)", palette = "YlGn", style = "cont") +
tm_shape(loss) +
tm_raster(title = "Forest cover loss 2001-2021",
style = "cat", breaks = c(0, 0.5, 1.5), palette = "red") +
tm_shape(loss_recent) +
tm_raster(title = "Forest cover loss 2022-2024",
style = "cat", breaks = c(0, 0.5, 1.5), palette = "purple") +
tm_shape(makay_union) +
tm_borders(col = "darkblue") +
tm_layout(legend.outside = TRUE) +
tm_add_legend(type = "line", labels = "Makay PA delimitation",
col = "darkblue") +
tm_shape(coordinates_sf) +
tm_dots(col = "black", size = 0.1, labels = coordinates_sf$location) +
tm_text(text = "location", ymod = -0.5) +
tm_view(set.view = c(10)) +
tm_scale_bar()Couvert forestier en 2000 et perte de couvert forestier depuis (source: GFC)